module basicmodule
	use myfuncs
	use dotops
	implicit none
	
	integer, parameter	:: nsim=10000

	logical				:: taumodel
	integer				:: tstind
	integer, parameter	:: tstNyb=0, tstxsi=1, tstqnt=2, tstPar=3, tstZpf=4 
	
	integer				:: k,n
	
	real, parameter		:: DGPxsimax=1.5	! max for size control grid
	real, parameter		:: qGEV=0.9

	
	real, allocatable	:: ys(:,:),tau(:),gammas(:,:,:)
	integer, allocatable	:: yks(:)
	real, allocatable	:: allstats(:,:,:)

	
	character(len=80)	:: dir="c:/dropbox/mystuff/TVextremes/out/"

	
!$omp threadprivate(ys,yks)	
	
	contains
	
	function DGPxsimin() result(val)
		real	:: val
		if(tstind==tstPar.or. tstind==tstZpf) then
			val=.03
		else
			val=-.5
		endif
	end function
	
	function getfilesuffix() result(val)
		character(:), allocatable	:: val
		if(.not. taumodel) then
			val="test"//convtos(tstind)//"_"//convtos(k)//"_"//convtos(n)
		else
			val="test"//convtos(tstind)//"_"//merge("tau","tvt",all(tau==1))//"_"//convtos(n)
		endif
	end function
	
	elemental function log1p(x) result(val)
		real, intent(in)	:: x
		real				:: val
		if(abs(x) > epsilon(1.0)**.33) then
			val=log(1+x)
		else
			val=x*(1-0.5*x)
		endif
	end function
	
	elemental function dlog1p(x) result(val)
		real, intent(in)	:: x
		real				:: val
		if(abs(x) > epsilon(1.0)**.33) then
			val=1/(1+x)
		else
			val=1-x
		endif
	end function
	
	elemental function expm1(x) result(val)
		real, intent(in)	:: x
		real				:: val
		if(abs(x) > epsilon(1.0)**.33) then
			val=exp(x)-1
		else
			val=x*(1+0.5*x)
		endif
	end function

	function getll_GEV(msx,Y,yk,ke) result(val)
        real    :: y(k),msx(3),yk,m,sig,xsi,val
		integer	:: ke
        real    :: zk
        m=msx(1); sig=msx(2); xsi=msx(3)
        zk=(yk-m)/sig
        if(1+xsi*zk<0) then
			if(ke==0 .and. xsi<0) then
				val=0
			else
				val=-1E50
			endif
            return
		endif
        if(ke>0 .and. 1+xsi*(y(1)-m)/sig<0) then
            val=-1E50
		else
			val=-ke*log(sig)-exp(-log1p(xsi*zk)/xsi)-(1+1/xsi)*sum(log1p(xsi*(y(1:ke)-m)/sig))
		endif
		if(.not. isfinite(val)) then
            val=-1E50
		endif
	end function
	
	function getdll_GEV(msx,Y,yk,ke) result(val)
        real    :: y(k),msx(3),yk,m,sig,xsi,val(3)
		integer	:: ke
		real	:: g,dgxsi,dldz
        real    :: zk,z,zxsi,dl1p
		integer	:: i
        m=msx(1); sig=msx(2); xsi=msx(3)
		zk=(yk-m)/sig
        if(1+xsi*zk<0) then
            val=0
            return
		endif
        if(ke>0 .and. 1+xsi*(y(1)-m)/sig<0) then
            val=0
			return
		endif
		val=0
		do i=1,ke
			z=(y(i)-m)/sig
			g=-log1p(xsi*z)/xsi
			dl1p=dlog1p(xsi*z)
			dgxsi=-g/xsi-(z/xsi)*dl1p
			dldz=-(1+xsi)*dl1p
			val(1)=val(1)-dldz/sig
			val(2)=val(2)-z*dldz/sig
			val(3)=val(3)+g+(1+xsi)*dgxsi
		enddo
		g=-log1p(xsi*zk)/xsi
		dl1p=dlog1p(xsi*zk)
		dldz=exp(g)*dl1p
		dgxsi=-g/xsi-(zk/xsi)*dl1p
		val(1)=val(1)-dldz/sig
		val(2)=val(2)-ke/sig-zk*dldz/sig
		val(3)=val(3)-exp(g)*dgxsi
		if(.not. all(isfinite(val))) then
!			print *,"problem in getdll_base"
            val=0
		endif
	end function
	
	function getconstll(msx) result(val)
		real	:: val,msx(3)
		integer	:: j
		val=0
		do j=1,n
			val=val+getll_GEV(msx,Ys(:,j),Ys(k,j),yks(j))
		enddo
	end function
	
	function getdconstll(msx) result(val)
		real	:: val(3),msx(3)
		integer	:: j
		val=0
		do j=1,n
			val=val+getdll_GEV(msx,Ys(:,j),Ys(k,j),yks(j))
		enddo
	end function

	function getscore(msx,y,ke) result(val)
		real	:: msx(3),y(k),val(3)
		integer	:: ke
		val=getdll_GEV(msx,Y,y(k),ke)
	end function
	
	subroutine setgammas
		integer	:: l,t,i
		real	:: eps(k)
		call rnset(13)
		do l=1,nsim
			do t=1,n
				call rnun(eps)
				eps=-log(eps)
				do i=2,k
					eps(i)=eps(i)+eps(i-1)
				enddo
				gammas(:,t,l)=eps
			enddo
		enddo
	end subroutine

	subroutine draw_Ys(msx,l) 
		real	:: msx(3),y(k),z(k)
		integer	:: l,t
		do t=1,n
			z=gammas(:,t,l)
			z=z**(-msx(3))
			z=(z-1)/msx(3)
			Y=msx(2)*z+msx(1)
			if(taumodel) then
				if(y(k)>tau(t)) print *,"warning: k too small in tau model"
				Y(k)=tau(t)
				yks(t)=count(Y(1:k-1)>tau(t))
			endif
			ys(:,t)=y
		enddo
	end subroutine
	
	function matinv(mat) result(val)
		use evcsf_int
		real	:: mat(:,:),val(size(mat,1),size(mat,1))
		real	:: evs(size(mat,1)),evecs(size(mat,1),size(mat,1))
		integer	:: j		
		call evcsf(mat,evs,evecs)
		do j=1,size(mat,1)
			val(:,j)=merge(1/evs(j),0.0,evs(j)>1E-7)*evecs(:,j)
		enddo
		val=matmul(val,transpose(evecs))
	end function
	
end module